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Abstract 

In this mini-review we summarize the progress of Lattice Boltzmann (LB) modeling and simulat- 
ing compressible flows in our group in recent years. Main contents include (i) Single-Relaxation- 
Time(SRT) LB model supplemented by additional viscosity, (ii) Multiple-Relaxation-Time(MRT) 
LB model, and (iii) LB study on hydrodynamic instabilities. The former two belong to improve- 
ments of physical modeling and the third belongs to simulation or application. The SRT-LB model 
supplemented by additional viscosity keeps the original framework of Lattice Bhatnagar-Gross- 
Krook (LBGK). So, it is easier and more convenient for previous SRT-LB users. The MRT-LB 
is a completely new framework for physical modeling. It significantly extends the range of LB 
applications. The cost is longer computational time. The developed SRT-LB and MRT-LB are 
complementary from the sides of convenience and applicability. 
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I. INTRODUCTION 



During the past two decades the lattice Boltzmann (LB) method has emerged as a com- 
Detitive scheme for simulating various nearly incompressible complex flows rangingfrom 



magnetohydrodynamics [2|, y j , to flows of suspensions |4(, flows with phase separation |5Hl4j. 
flows through porous media flQ, etc. With increasing the Mach number, the compress- 
ibility of flow becomes more pronounced. Such high speed compressible flows are ubiquitous 
in explosion physics, aerophysics and astrophysics, etc. Up to now, the LB modeling and 
simulating of compressible flows, especially those with shocks and/or discontinuities, is still 
a challenging issue. 

Given the great importance of shocking and detonation in many fields of physics and 



engineering constructing LB models for high speed compressible flows has been 

attempted since the early days of LB research To proceed, we first discuss the most 
fundamental problem "what is LB ?". The views are not exactly the same in papers by 
different authors. Since having different knowledge backgrounds and working in different 
fields, different authors may use LB to solve different problems and focus on different sides 
of LB. Understandably, even for the same author, the views will be updated with extending 
research experience. Globally speaking, the views on LB can be classified into two categories. 
The first category regards LB as a new scheme for simulating hydrodynamic equations such 
as the Euler equations and Navier-Stokes equations. The second category regards LB as a 
kind of new model of physical systems. Physical model construction and numerical method 
design are the first two steps for numerical study on any physical problems. Compared with 
numerical methods, the physical model construction is the first step and more fundamental. 
Only after the physical model is fixed can the corresponding numerical method be estab- 
lished. Clearly, the first kind of view starts LB research from the second step, numerical 
method design. It does not consider the improvements of the physical modeling. In other 
words, it assumes that the original hydrodynamic equations are sufficiently exact for model- 
ing the problem under consideration. The second kind of view puts LB research on the more 
fundamental step, physical modeling. For this view the numerical method is the second 
important issue. It accepts any reasonable numerical methods no matter they are new or 
traditional. The second kind of view aims at physical problems. The point of the second view 
is that, compared with the traditional hydrodynamic equations, the LB framework contains 
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more physical components. The theoretical reasons are as below. The LB model is based on 
the Boltzmann equation which is one of the most fundamental equations in non-equilibrium 
statistical physics. It naturally inherits some intrinsic characteristics of the latter. Accord- 
ing to the Chapman-Enskog analysis, one can expand the distribution function around its 
equilibrium as Taylor series in the Knudsen number. When the Knudsen number approaches 
zero, the system is nearly in equilibrium state, the deviation from equilibrium is negligible, 
the LB model corresponds to or recovers the Euler equations. When the first order terms in 
Knudsen number have to be accounted and the second order terms are negligible, in other 
words, when the system slightly deviates from the equilibrium, the LB model corresponds to 
or recovers the Navier-Stokes equations. When the system deviates more from equilibrium 
and the second order terms in Knudsen number have to be taken into account, the LB model 
is beyond the Navier-Stokes description. The theoretical framework of LB is self-adaptive 
for describing complex systems where the deviations from equilibrium are spatially and tem- 
porally varying. From the view of modeling precision on detailed dynamics, it is less than 
Molecular Dynamics(MD). It adopts the concept of distribution function. It is generally 
considered as a kind of mesoscopic modeling. For continuum system, the LB should give 
the same results as those of hydrodynamics equations. For non-continuum systems such 
as the boundary layers where the Knudsen number is high, the LB should give the same 
results as those of other mature methods such as MD or Monte Carlo(MC). In between the 
two kinds of limiting cases, the hydrodynamic equations are not valid, the MD and MC are 
reasonable but not practical due to the huge quantity of computations. For such cases, the 
LB modeling and simulation still work. Its results should be checked by physical principles 
and analyses. Just as in traditional Computational Fluid Dynamics(CFD) where different 
discretization schemes work for different problems, for different systems one should compose 
or choose different LB models. 



In 1992 Alexander et al [20|] proposed a compressible LB model where the main skill is to 
introduce a flexible sound speed so that the Mach number may become higher. This model 
works only for nearly isothermal compressible systems. In 1999 Yan et al 2lJ proposed a LB 
scheme for compressible Euler equations. In this model a Discrete Velocity Model(DVM) 



with three energy levels is used. Sun et al 22|, |23[ proposed an adaptive LB model where the 



particle velocities vary with the Mach number and internal energy. The model partly frees 
the particle velocity from fixed values. It works for more extensive systems compared with 
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previous LB versions. Its two-dimensional and three-dimensional versions were published 
in 1998 and 2003, respectively. The evolutions of all those models follow the traditional 
"propagation + collision" mode. All of them belong to the standard LB models. Due to 
the inconvenience of application and/or numerical instability problems, few physical results 
based on those models can be found. 

For modeling and simulating compressible flows, an alternative way is to use the Finite- 
Difference (FD)-LB method. Tsutahara group 24j-|28| in Kobe university proposed several 
FD-LB models in recent years. The FD-LB model frees the combination of spatial and 
temporal discretizations. The sizes of particle velocities are flexible. So it is much more 
convenient to meet the requirements for simulating compressible fluids. The FD-LB scheme 



was then extended to the case of binary fluids [29|, [30|. But numerical instability problem 
blocks its practical applications to systems with a Mach number being larger than 1. In fact, 
as for the numerical instability problem, many attempts have been made. Typical examples 
are referred to the entropy LB model 3jJ,|32|, FIX-UP scheme 33|], flux-limters approach [341] . 
etc. But most of the discussions were still focused on systems with small Mach numbers. 

To model and simulate high speed compressible flows, especially those with shocks, our 
group developed two schemes in recent years. The first is to introduce additional viscosity 



and improve the discretization of spatial and temporal derivatives [35|-|40|. This scheme 
does not change the framework of the original LB model. The second is to develop Mul- 
tiple Relaxation Time(MRT) LB models 4ll445|. The framework is changed in the second 
scheme. The first scheme is based on the following facts, (i) The numerical fluid particles 
do not distinguish the original viscosity and additional viscosity, (ii) Introducing additional 
viscosity is equivalent to modifying the relaxation time from some sense, (iii) Better tem- 
plate of discretization may damp the numerical anisotropy. Our improved models work for 
both high speed and low speed flows. So, they make it possible to simulate stable shocks in 
compressible fluids. The first scheme is based on the original Bhatnagar-Gross-Krook(BGK) 
model. It is a remedy under the original framework. 

The rest of the paper is structured as follows. We first introduce a few improved LB 
models based on the first scheme in section II. The MRT scheme is reviewed in section 
III. Section IV shows two typical applications, LB studies on Richtmyer-Meshkov(RM) and 
Kelvin- Helmhotz(KH) instabilities. Section V summarizes the present paper. 
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II. SRT MODEL SUPPLEMENTED BY ADDITIONAL VISCOSITY 



Among the two-dimensional FD-LB models for compressible flows, the one by Kataoka 



and Tsutahara 



24( is typical. It has very simple and strict theoretical background, uses 
a DVM with only 9 components. The specific heat ratio is flexible. But the numerical 
instability blocks its application in supersonic flows. Therefore, our first LB model for high 
speed compressible flows is created by improving the Kataoka-Tsutahara(KT) model. 
The LB kinetic equation with BGK approximation reads, 

% + v ia ^ = -[fi eq -M (i) 

at ox a t 

where /, {fi q ) is the discrete (equilibrium) distribution function; v, is the i-th discrete 
velocity, % — 0, • • • , N — 1; N is the total number of the discrete velocity; index a = 1, 2, 
3 corresponding to x, y, and z, respectively; r is the relaxation time determining the speed 
of approaching equilibrium. Sometimes, r is rewritten as er', where e is a dimensionless 
number, the Knudsen number. The original KT model corresponds to the complete Euler 
equations 



dp d(pu a ) = Q 
dt dx a 

d(pug) + d(pu a u p ) + d^ = Q , 2 v 
dt dxp dx a 

Ft {E + \ pul) + lk [Ua{E + \^ + P)] = °' 

when the knudsen number e approaching zero. Here p, u, P (= pT), E(= pT/fy — 1)) are 
the hydrodynamic density, flow velocity, pressure and internal energy, respectively; T is the 
temperature and 7 is the specific-heat ratio. To make 7 flexible, a constant, b = 2/(7 — 1), 
is introduced. The following constraints are needed for this model, 

N-l N-l 

p=T,fi q = T,f» ( 3 ) 

i=0 i=Q 
N-l N-l 



p Ua = f? q V ia = Y f iVia ' ( 4 ) 

i=0 i=0 
N-l N-l 

p(bRT + ul) = J2 f!>l + ^) = E M< + ^ ( 5 ) 



i=0 i=0 



FIG. 1: Schematic figure of the discrete velocity model. 



N-l 



i=0 
AT-1 



P 



[(b + 2)RT + 4] u a = Y,ft>% + ^ 



i=0 



where rji is another variable introduced to make specific-heat ratio flexible. 
In the two-dimensional case, the KT DVM has nine components. It reads 



(0,0), 



Ci cos 



• 7r(i+l) - 



sin 



{va,Vi2} 



■7r(i+l) - 



1,2,3,4 



c 2 [cos vr(^i + i), sinvr(^l + J)], % = 5, 6, 7, 



(6) 
(7) 



(8) 



(9) 



r) , i = 

0, f = 1,2,...,8 ' 

A schematic figure of the distribution of the discrete velocities is shown in Fig.l, where 
c\ and C2 are constants which should not depart faraway from the flow velocity u. C2 is 
generally chosen 1.0 ~ 3.0 times of c\. 

The local equilibrium distribution function is computed by 



f- 9 = p(Ai + BiV ia U a + DiUaViaUpVip), % = 0,1, • • • ,8, 



(10) 



where 
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(11) 
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Parameters r]o, c\ and C2 are independent in this DVM. 770 influences fl q via the expansion 
coefficient A4. In the original KT model, the usual FD scheme with first-order forward in 
time and second-order upwinding in space is used. 

To make practical the LB simulation to the supersonic flows, we propose an alternative FD 
scheme combined with an additional dissipation term to overcome the numerical instability 
problem. The LB equation ([1]) can be regarded as non-dimensional. In this work, we consider 
t = er' and set the time step At to be numerically equal to the Knudsen number e. Thus, 
from Eq.flT]) we have 

/,(x,t + At) - /,(x,t) + v ia ^p^At = - [/^(x,t) - ffat)] . (13) 

ox a r 

In Eq. (fT3j) r' has been written as r for simplicity. The spatial derivative dfi/dx can be 
calculated by 

„ ^ n dfi /3fj(x + Ax, t) + (1 - 2(3) f(x, t)-(l- P)fi(x ~ Ax, t) 
If ^"°' & = Ax ; (14) 

„ ^ n dfj _ (1 - m{x + Ax, t) - (1 - 2(3) f(x, t) - 0Mx - Ax, t) 
l±«fa<0, ^ . (15) 

In Eqs. p4p and f|T5|) . < (3 < 0.5. If /3 takes zero, then they are no other than the first-order 
upwind scheme in space; if /3 takes 0.5, they recover to the general central difference scheme. 
df/dy can be calculated in a similar way. Actually, Eqs. ffT4l) and (fl5|) can be rewritten as 

« ^ n d fi fi(x,t) - f(x- Ax,t) 

If V - ~ °' = Ax (16) 

(3Ax[fj(x + Ax, t) + fj(x- Ax, t) - 2fj(x, t)} 

+ Ax 2 
Tf ^ n dfi f(x + Ax,t) - fi(x,t) 

lf^<0, ^ = ^ (17) 

/3Ax[/i(x + Ax, t) +fi{x- Ax, t) - 2f(x, t)} 



Ax 2 

The second terms in the Right-Hand-Side(RHS) of Eqs. flTE]) and ( ITT)) can be regarded as 
some kind of additional viscosities which can reduce some unphysical phenomena such as 



wall-heating, but they are not enough. Additional dissipation term is needed. The final LB 
equation reads 



/«(x, t + At)- /,(x, t) + v ia ^^-At - \ 92 f*2 t) At = ^^(x, t) - /,(x, t)} (18) 

° 0=1 a 

where A« is a small number not varying in space or time. The second-order derivative 9 q^'^ 
can be calculated by the central difference scheme. In our simulations Ax = Ay and the 
parameter /3 is generally chosen to be 0.25 if not particularly claimed. How to choose the Aj 
is the key problem. Analysis by the software, Mathematica, and numerical tests show that 
we can choose Aj around the following way, 

c\Ax, i = 

ci Ax/10, i = 1,2,3,4 • (19) 
0, 2 = 5,6,7,8 

The improved model is validated by well-known benchmark tests. Simulations on Rie- 
mann problems with very high ratios (1000:1) of pressure and density also show good ac- 
curacy and stability. Regular and double Mach shock reflections are successfully simulated. 
It should be commented that, since using constraint, At = e, such a model can only be 
regarded as a new scheme to simulate the Euler equations. The added viscosity terms can 
be regarded as a kind of slight remedy to the traditional hydrodynamic model. 



In 2008 Gan, Xu, Zhang, et al [36j developed a LB model for high speed compressible 
flows. In this model, the constraint, At = e, is eliminated. Therefore, it can be regarded 
as a mesoscopic new model. In the continuum limit it corresponds to the Navier-Stokes 
equations. The model is composed of three components: (i) the DVM by Watari and 
Tsutahara 2fJ, (ii) a modified Lax-Wendroff FD scheme where reasonable dissipation and 
dispersion are naturally included, (iii) additional viscosity. The improved model is convenient 
to compromise the high accuracy and stability. The included dispersion term can effectively 
reduce the numerical oscillation at discontinuity. Shock tubes and shock reflections are 
used to validate the new scheme. In our numerical tests the Mach numbers are successfully 
increased up to 20 or higher. In Fig. 2 we show a simulation result on double Mach reflection 
by the improved model. The initial pressure ratio here is high. A planar shock is incident 
towards an oblique surface with a 30° angle to the direction of propagation of the shock. A 
uniform mesh size of 500 x 200 is used for the numerical simulation. The conditions for both 
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FIG. 2: Contours of density (a), temperature (b), and u x (c) of the double Mach reflection problem 
at the time t = 7.5 x 10 -3 . The units of the x- and y- axes are both 0.001. 



sides are: 



' 400 
^ 67 



(P; 'U'xi U"yi 2"") |x,j/,0 

,13.3cos30°, - 13.3 sin 30°,89.2775), if y > h(x,0) 



(20) 



(2.0,0.0,0.0,0.5), if y < h{x,0) 

where h(x, t) = a/3(x — 80Aa:) — AOt. The reflecting wall lines along the bottom of the 
problem domain, beginning at x = 0.08. The shock makes a 60° angle with the x axis and 
extends to the top of the problem domain at y = 0.2. At the top boundary, the physical 
quantities are assigned the same values as on the left side for x < g(t) and are assigned the 
same values as on the right side for x > g(t), where g(t) = 80 Ax + ^3(0.2 + 40t). The 
computed density, temperature and flow velocity along the x-direction are shown in Fig. 2, 
where complex characteristics, such as oblique shocks and triple points, are well captured. 
In this model the ratio of specific heat is fixed on an unphysical constant 2. Later, 



Gan, Xu, Zhang, et al studied a model for flexible specific heat ratio 



371 ] . For higher 



computational efficiency, Chen, Xu, Zhang, et al proposed a model where the number of 
discrete velocity decreases from 65 to 16 [38:] . They simulated the reaction of shock wave 
on a bubble or ball, etc. In 2010 they present a three-dimensional LB model for high Mach 
number compressible flows. Figures 3(a) and 3(b) show our successful LB simulations of 
shock wave reactions on bubble and on ball, respectively, where only the density isosurfaces 
are shown. In both Figs, (a) and (b), the upper plot shows the initial state, and the lower 
one shows a snapshot in the shocking procedure. The added additional viscosity makes the 
scheme more consistent with the physical system and more convenient to satisfy the von 
Neumann stability condition. Among the discussions on LB model with additional viscosity, 



the application of flux limiters is also investigated 40] . In the reference with flux limiters 



40J Gan, Xu, Zhang, et al also introduced an improved BGK model to break the fixed- 
Prandtl- number barrier. It is meaningful to briefly review the scheme for this improvement. 

In the SRT model, both the viscosity and heat conductivity coefficients are proportional 
to the relaxation time r. As a result, the Pr is fixed to 

Pr= ^ = l. (21) 

K 

The control of Pr may be achieved by modifying the BGK collision term as below: 

+ v« ■ = ~\ [hi - (1 + Ar)#] , (22) 
where A takes the following form 

A = A + B(y kl - u) 2 . (23) 

Contributions of the new term Af^ in Eq. (|22|) to the mass, momentum, and energy equations 
are 

J2 A fk1 = (A + 2BT)p = 0, (24) 

ki 

A fki^ a = (A + 2BT)pu a = 0, (25) 

ki 

E \ A fki v l = P(A + 2BT)(T + ^) + 2pT*B = 2pT 2 R (26) 

ki 

We require that Eq.( l22l) recovers the Navier-Stokes equations in the following form, 

dp d{pu a ) 



dt dr r 



(27) 
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d(pu a ) d(pu a U/3 + P5 a p) d r ,duf3 du a , ->\^ u i 



+ _ = 

at or/3 orp or a or^ c/r 7 

where \i = pTr is the viscosity, k is the heat conductivity, k is required to be k = c p pT(r+q), 
where c p = jCy = 7/(7 — 1) is the specific-heat at constant pressure. It is clear that a new 
coefficient q is introduced to make the Prandtl number flexible. By using Eqs. (l24l) - (l26j) it is 
easy to find coefficients in Eq. fl23|) with the following form 



A = -2BT, B = -±—d a {c p pTqd a T}. (30) 
2pl z 

Therefore, the modified BGK collision term changes the heat conductibility in the energy 
equation from k = c p pTr to k — c p pT(r + q). Consequently, the Prandtl number is changed 
to 

Pr = (31) 

r + q v 1 

Figure 4 shows a validation example of such a scheme for flexible Prandtl numbers based 
on the SRT model. The figure shows the comparison of LB results with theoretical solutions 
for thermal Couette Flows. Fig. (a) is for the temperature profiles in steady state for various 
Prandtl numbers. Fig.(b) shows the velocity profiles for Pr = 5.0 at various times. For 



more details the readers can refer to Ref. 



401 ] . Such a scheme makes a significant remedy 



from the side of physical modeling. It is easy to find that such a scheme can also be used to 
change other transport coefficients such as the viscosity. It is also meaningful to mention that 
among the moment relations required by each LB model, only for the three, the definitions 
of density, momentum and energy, the equilibrium distribution function f^ q can be replaced 
by the distribution function If we replace f^ q by fi in RHS of any other required moment 
relations, the value of RHS will have a deviation from that of the left hand side. This 
deviation may work as a measure for the deviation of system from its equilibrium. For 
example, the following Ai 

JV-l N-l 

A l = fi V i<* V % ~ f? V icV%, (32) 
j=0 i=0 

presents a measure for how much the system deviates from its equilibrium for cases without 
using the constraint At = e. 
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(a) (b) 

FIG. 3: Density isosurfaces of shocked bubble (a) and shocked ball (b). In (a) or (b) the upper 
plot shows the initial state, the bottom one shows the density configuration during the shocking 
procedure. 
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FIG. 4: Comparison of LB results and theoretical solutions for thermal Couette Flows, (a) 
Temperature profiles in steady state for various Prandtl number, (b) Velocity profiles for Pr = 5.0 
at various times. 

III. MRT MODEL 

It is known that different motion modes generally approach their equilibria in different 
velocities. But in the SRT BGK model, the speeds of all discrete distribution functions 
approaching the equilibria are determined by a single relaxation time r. That means r is 
an averaged relaxation time of all kinds of motion modes. The best merit of this treatment 
is that it is simple and keeps the most fundamental conservation laws. This BGK model 
has been successfully applied in various fields. But with increasing the Mach number and 
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Reynolds number, the problem of numerical instability becomes more serious. At the same 
time, the Prandtl number effect is a key issue in many ffuidic systems. Facing with all these 
requirements and challenges, people began to reevaluate this simple averaging treatment. 

The numerical instability of LB simulation is still a difficult problem nowadays. Roughly 
speaking, the possible reasons come from two sides, the physical modeling and the dis- 
cretization scheme. It has been indicated that untying the motion modes which should be 



independent is helpful for improving the numerical stability 



et al 



48 



46to0| . Succi, et al [461 ]. Luo 



- 50| | and many others have made significant contributions in constructing MRT LB 
models. Those MRT models are mainly within the framework of the standard LB model and 
work for isothermal systems with low Mach number. In recent years our group proposed two 
schemes to compose MRT model for high speed compressible flows. These schemes are for 
the framework of the FD-LB model. The finished works focus still on the two-dimensional 
cases. 

In the MRT LB formulation, the collision step is first calculated in the kinetic moment 
space spanned by a suitable set of N kinetic moments of the distribution function /j. Then, 
the propagation step is performed back in the discrete velocity space spanned by the N 
discrete velocities Vj. In contrast to the SRT model, the MRT version caters for more 
adjustable parameters and degrees of freedom. The relaxation rates of the various kinetic 
moments due to particle collisions may be adjusted independently. The MRT LB equation 



has the following form, 



dfi dfi 

-gj- + Via g^r = _S * fc lf k ~ fh\ ' ( 33 ) 



where S is the collision matrix. The equation reduces to the usual lattice BGK equation if 
all the relaxation parameters are set to be a single relaxation time r, namely S = -I, where 
I is the identity matrix. The discrete distribution functions /j and f^ q can be rewritten as 
the following matrixes: 

f =(/i,/ 2 ,--- Jn) T , (34a) 
f e " = (/f , f?, ■ ■ • , f N q ) T , (34b) 

where T is the transpose operator. Given a set of discrete velocities v$ and corresponding 
distribution functions /j, we can get a velocity space S spanned by discrete velocities Vj 
and a moment space S M spanned by moments of particle distribution function fa. The 
moments of particle distribution function reads f = ^/i,/2, •" > f^j , where /$ = m^fj, 

13 



rriij is an element of the matrix M and is a polynomial of discrete velocities. Obviously, 
the moments are simply linear combination of distribution functions fi, and the mapping 
between moment space and velocity space is defined by the linear transformation M, i.e., 
f = Mf, f = M _1 f , where M = (m 1; m 2 , • ■ ■ ,m^) T , m, = m i2 , • ■ ■ , rrijjv)- 

Since the collision step is first calculated in the moment space and then mapped back to 
the velocity space. So, the MRT LB equation can be described as 

^ + v ia ^ = -M u % k (f k -f k q ), (35) 

where S = MSM -1 = diag(si, s 2 , ■ ■ ■ , sjv) is a diagonal relaxation matrix. f e { q is the 
equilibrium value of the moment /j. The moments can be divided into two groups. The first 
group consists of the moments locally conserved in the collision process, i.e. fi = f^ q . The 
second group consists of the moments not conserved, i.e. fi ^ f^ q . The equilibrium f^ q is 
a function of conserved moments. It is clear that the first group includes the density, the 
momentum and the energy. 



A. MRT model based on group representation theory 



Now we briefly review the first MRT LB model proposed in our group l4l|. Our first 
MRT model is developed from the SRT version by Kataoka and Tsutahara [25]. The DVM 
can be expressed as: 



eye: (±1,0), forl<i<4, 
eye : (±6, 0) , for 5 < i < 8, 
V2(±l,±l) , for 9 < z < 12, 



(36) 



k ^(±1,±1), for 13<z< 16, 
where eye indicates the cyclic permutation, (see Fig. 5) 
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FIG. 5: Distribution of Vj Q for the discrete velocity model. 
1. Construction of transformation matrix M 

The transformation matrix M is constructed according to the irreducible representations 
of SO (2) group: 

1- 

cos 9, sin 9, 

sin 2 9 + cos 2 9, cos 29, sin 29, 

cos6 l (sin 2 9 + cos 2 9), sin ^(sin 2 9 + cos 2 9), cos 3$, sin 36*, 
(sin 2 9 + cos 2 9) 2 , cos 40, cos 2fl(sin 2 9 + cos 2 9), 
sin 2#(sin 2 # + cos 2 #), 

cos3#(sin 2 9 + cos 2 9), sin3#(sin 2 9 + cos 2 9), ■ ■ ■ 

Let Vi X and Vi y play the roles of cos 9 and sin#, respectively. Then we define mu = 1, 
m 2i = v ix , m 3i = v iy , m 4i = (v 2 x + v 2 y )/2, m 5i = v 2 x -v 2 y , m 6i = v ix v iy , m 7i = v ix (v 2 x + v 2 y )/2, 
mgi = v iy [v 2 ix + v 2 y )/2, m 9i = v ix {v 2 x - 3v 2 y ), m Wi = v iy (3vf x - v 2 y ), m m = (vf x + v 2 y ) 2 /4, 
m 12i = vf x - Qv 2 x v 2 y + vf y , mis, = (i& + v 2 y )(v 2 x - v 2 y ), m Ui = {v 2 x + v 2 iy )v ix v iy , m 15i = 
v ix(v 2 ix - 3v 2 y )(v 2 x + vfy), m Wi = v iy (3v 2 x - v 2 iy ){v 2 ix + v 2 y ), where i = 1, ■ ■ ■ , 16. 

For two-dimensional compressible models, we have four conserved moments, density f\ = 
P = Y.fi m Ui momenta f 2 = j x = pu x = ]T f i m 2i and f 3 = j y = pu y = Yl fi m 3i, and 
energy = e = ^ jivn^. To be consistent with the idiomatic expression of energy, in the 
definitions of m^, mji and mgj, a coefficient 1/2 is used. Similarly, a coefficient 1/4 is used 
in the definition of m Ui . The components of transformation matrix M are shown in table I. 
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TABLE I: Transformation matrix of MRT-LB for compressible fluids. 



1 


mu 


m 2 i 


m 3i 




m 5i 


m 6 i 


777,7; 


m 8 i 


mgi 


mioi 


mm 


m-42i 


mi3i 


777,14, 


m 15i 


mm 


1 


1 


1 





i 

2 


1 





1 

2 





1 





1 

■1 


1 


1 





1 





2 


1 





1 


1 

2 


-1 








1 

2 





-1 


1 

1 


1 


-1 








-1 


3 


1 


-1 





1 

2 


1 





1 

2 





-1 





1 

4 


1 


1 





-1 





4 


1 





-1 


1 

2 


-1 








1 

2 





1 


1 

4 


1 


-1 








1 


5 


1 


6 





IS 


36 





108 





216 





324 


1296 


1296 





7776 





6 


1 





G 


18 


-36 








108 





-216 


324 


1296 


-1296 








-7776 


7 


1 


-6 





18 


36 





-108 





-216 





324 


1296 


1296 





-7776 





8 


1 





-6 


18 


-36 








-108 





216 


324 


1296 


-1296 








7776 


9 


1 


V2 


V2 


2 





2 


2V2 


2^2 


-4\/2 


4V2 


4 


-16 





8 


-16\/2 


16^2 


10 


1 


-V2 


V2 


2 





-2 


-2^2 


2V2 


4^/2 


4V2 


4 


-16 





-8 


16V2 


16V2 


11 


1 


-V2 


-V2 


2 





2 


-2^2 


-2V2 


4^/2 


-4V2 


4 


-16 





8 


16^/2 


-16\/2 


12 


1 


V2 


-V2 


2 





-2 


2^2 


-2^2 


-4^/2 


-4V2 


4 


-16 





-8 


-16^2 


-16\/2 


13 


1 


3 

V2 


3 

V2 


9 
2 





9 
2 


27 


27 
2V2 


27 


27 
V2 


si 
4 


-81 





si 
2 


243 
V2 


243 
V2 


14 


1 


3 

V2 


3 

V2 


9 
2 





9 
2 


27 


27 
2^2 


27 
V2 


27 
V2 


SI 
4 


-81 





84 
2 


243 
V2 


243 
V2 
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1 


3 

V2 


3 

V2 


9 

2 





9 
2 


27 
2\/2 


27 
2V2 


27 
V2 


27 
V2 


SI 
4 


-81 





SI 
2 


243 
V2 


243 
V2 


16 


1 


3 

V2 


3 

V2 


9 
2 





9 
2 


27 
2V2 


27 
2V2 


27 
V2 


27 
V2 


SI 
4 


-81 





84 
2 


243 
V2 


243 
V2 



2. Determination of f^ q 

The second group components of f^ q are chosen in such a way that in the continuum 
limit the MRT LB model recovers the Navier-Stokes equations. To that end, we perform the 
Chapman-Enskog expansion on the two sides of Eq. (l33|) . We use the following multiscale 
expansions: 

f i = fi 0) +fP + f?\ (37a) 
Odd . „ , 

st = ST + (37b) 

dx dxi ' ^ 
where is the zeroth order, ff~\ ^-and are the first order, f- 2 ^ and J^- are the second 
order terms of the Knudsen number e. Equating the coefficients of the zeroth, the first, and 



1(3 



the second order terms in e gives 



C = n\ (38a 



1 dt\ dx 



(0) _ f eq 
J i i 



d d , f(0) _ (i) 



la 



In the moment space they are 



f?> = ST, (39a) 

^ + ^£^ = -^ (39b) 
^' + (^ + = -S„/f, (390) 

where E Q = M(f ia I)M _1 . 

The equilibria of the moments in the moment space read f eq = 
(p,j x ,jy,G,f^ 9 ,fQ q ,--- , fil) T ■ The first and second order deviations from equilibria are 
defined as : f« = (0, 0, 0, 0, / 5 (1) , / 6 (1) , ■ • • , f^f and f< 2 ) = (0, 0, 0, 0, jf, ff\ ■ • ■ , f^f , 
respectively. Via some algebraic treatments, we obtain 

?E + -k. + -hi = o f40a) 
dt dx dy 

If choose j? = Ul ~ f,)IP, ft' = J,h/P, ff = (e + /f = (e + P)jjp, !? = 

(A ~ 3f w ) 3 Jp\ US = - f„)3jp 2 , fS = 2e 2 /p - (£ + j}) W, /S = (6pe - 2f - 



X 



- 3 2 y )/p\ ffi = (6pe - 2f x - 2f y ) 3xh /p\ fS = fit = fit = 0, the MRT LB model 
recovers the following Navier-Stokes equations: 



?E + ?h. + ?h. = o (41a) 

dt dx dy 

dj x i 9 ( .2 1 _\ ® (■ ■ i \ ®P , ^ r ^ m s/m , ^ r 



<9t dx dy dx dx dx dy dy dx dy 

§ + #■ (Wp) + #■ (£/,) = + 1- + £)] - £w£ - £ 

ot ax dy x y J dy dx dx dy dy dx dy 



^)]> ( 41b ) 
)], (41c) 
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dx 1 dx 2 y dx x dx x dy y dy 

d dT A 2 , du x du x du y du y 
+ dy- [X ^ + ~2 ~dy ~ ^ + ^ + ^ )] ' (4W) 

where /i s = pRT/s 5 , p v = pRT/s e , Ai = 2pRT/s 7 , A 2 = 2pRT/s$. It is noted that the 
definitions of /ff , /ff, have no effect on the recovered macroscopic equations. When 
p s = p v = p, Ai = A2 = A, the above Navier-Stokes equations reduce to 

<9t <9x Q <9x a dx x af3 

In Fig.6(a) we show an example of stability comparison for the new MRT model and its 
SRT version. The abscissa is for kdx, and the vertical axis is for |o;| max which is the largest 
eigenvalue of coefficient matrix Gij. The macroscopic values in stability analysis are chosen 
as follows: (p, u x , u y , T) = (2.0, 10.0, 0.0, 2.0). The relaxation time in SRT is r = 10 -5 , while 
the collision parameters in MRT are s 5 = 6500, s 7 = s$ = 9 x 10 4 , s 9 = 8 x 10 4 , S\ 3 = 7 x 10 4 , 
S14 = 8 x 10 3 , S15 = 2.5 x 10 4 , the others are 10 5 . In this case, the MRT scheme is stable, 
while the SRT version is not. It is clear that, by choosing appropriate collision parameters, 
the stability of MRT can be much better than the SRT. 

Figure 6(b) shows the comparison of MRT LB results and exact ones for the well-known 
Colella explosion wave problem. For the problem, the initial condition is 

(p,u x ,u y ,T)\ L = (1.0,0.0,0.0,1000.0), x<0. 
(p,u x ,u y ,T)\ R = (1.0,0.0,0.0,0.01), x > 0. 

This is a strong temperature discontinuity problem that can be used to study the robust- 
ness and precision of numerical methods. Figure 6(b) gives density, pressure, velocity and 
temperature results at t = 0.1. Symbols are for simulation results. Here, the parameters 
are s 7 = s 8 = 5 x 10 4 , s n = s 13 = 5 x 10 5 , other values of s still adopt 10 5 . The success 
of the simulation shows that the MRT model is applicable to simulate strong temperature 
discontinuity problem, and confirms the robustness and precision of the model. 
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FIG. 6: (a) Stability comparison for the new MRT model and its SRT version, (b) The MRT 
simulation results and exact solutions for the Colella explosion wave at time t = 0.1. 

Two points should be commented here. The first is that the better stability is not the only 
or most important advantage of MRT over SRT. From the view of physical modeling, the 
SRT is only a special case of the MRT. The second is that the above MRT LB model works 
well for shocked compressible fluids where the shocking procedure is much faster than the 
transportation processes. To work also well for more general cases, the collision operators 
of the moments related to the energy flux should be modified as below 45] , 



S 77 (/ 7 - ft) => S 77 (/ 7 - ft) + (s 7 /s 5 - l)pTu x (^ - ^) + (s 7 /s G - l)pTu y (^ + ^), 

y x 

S 88 (/ 8 - ft) =► S 88 (/ 8 - ft) + (s 8 /s 6 - l) P Tu x (^L + _ {S8 / S5 _ l)pTu y {^ - 



dx dy 



dx dy 
(45) 



After the modification the coefficients of viscosity in energy equation (141d|) are consistent 
with those in momentum equations ( 141b[) -( l4Tcl) . 



B. MRT model based on moment relations 



In the original KT model, besides Eqs. (E])-([7]), the local equilibrium distribution function 
f^ q is required to satisfy the following two additional moment relations: 



p [RT (u a 5p x + up5 ax + u x 5 a p) + u a upu x ] = ^ ft q Vi a Vipv ix , 



(46a) 



19 



p { (b + 2) R 2 T 2 5 a p + [(b + 4) u a up + u\5 a p\ RT + u 2 x u a up) = /f 9 (u? x + 77?) u ia v i/3 

(46b) 

The local equilibrium distribution function f? q is calculated via the following polynomial: 



fi 9 = Pl a oi + a u T + a 2iT + (a 3i + a 4i T) w a + a 5i u a u^ 

+ (6oi + + hiU 2 a ) upVip + (d i + d Xi T + d 2 i?4) upv il3 u x Vi 
+eiU a v ia UfsVii3U x Vi X }, (47) 

which is of the flow velocity up to the third order. The coefficients aoi , . . ^ej (i — 1, ... ,16 
) in the distribution function f^ q are referred to the original publication 25]. 



1 . Construction of transformation matrix M 



In this MRT model the moments are chosen according to the seven required moment 



relations 
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43| . The RHS of the seven equations indicate seven monomials: 1, v ia , v 2 a + r] 2 



ViaVi/3, {vfp + r] 2 )v ia , Vi a VipVi X , {y 2 x + r] 2 )v ia Vi/3. Three possibilities arise from the monomial 
ViaVip: (a) a = /3 = x, v ia v^ = v 2 x , (b) a = (3 = y, v ia v if} = v 2 y , (c) a = x, (3 = y, 
ViaVip = v ix v iy . "(a) + (b)" gives (v 2 x + v 2 y ), "(a)-(b)" gives (v 2 x - v 2 y ). Through such a simple 
combination of these monomials, we can compose the transformation matrix M as below: 
m u = 1, m 2i = v ix , i 

m 7i = v ix v iyi m 8i = "ix\"i x ' "i y ' <l%)i ,,U W ~ "iy\"ix 1 "iy 1 'HIl " D 1UJ — "ix\"i x ' "iyji 

mm = v iy (v 2 x +v 2 y ), mvH = v ix (v 2 x -v 2 y ), m m = v iy (v 2 x -v 2 y ), m 1M = (v 2 x +v 2 y )(v 2 x +v 2 y +r]f), 
m 15i = v ix v iy {v 2 x + v 2 iy + 77?), m m = {v 2 x - v 2 y ){v 2 x + v 2 y + r] 2 ), where % = 1, ■ • • , 16. The 
components of transformation matrix M are shown in table II. 



7;? — 7;? 

w ix iy> 



ixi " u 6i u iyi ■■•"a ~tx ' "iy ' ' / % 1 ""oi ~ix ' "ly 

vU v l + v l + Vf), m 9i = v iy (v 2 x + v 2 + 77?), m Wi = v ix (vf x + v?) 



It should be pointed out that, different from the other MRT models for isothermal fluids, 
the transformation matrix M should not be based upon a Gram-Schmidt orthogonalization 
procedure. 

2. Determination of f^ q 

The procedure of determining f eq is similar to that for the first MRT LB model in this 
paper. But the results are significantly different. Our choice for this model is as below: 
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TABLE II: Transformation matrix of MRT-LB model for compressible flows with flexible specific- 



heat ratio. 



i 


mu 


m 2 j 


m 3i 




m 5i 


m 6i 


m 7i 


m 8i 


m 9i 


mm 


mm 


m\2i 


m 13i 


mui 


m\u 


mm 


1 


1 


1 





29 
4 


1 


1 





29 
4 





1 





1 





29 
4 





29 
4 


2 


1 





1 


29 


1 


-1 








29 





1 





-1 


29 
4 





29 
4 


3 


1 


-1 





29 
4 


1 


1 





_29 
4 





-1 





-1 





29 
4 





29 
4 


4 


1 





-1 


29 

1 


1 


-1 








_29 

A 
1 





-1 


1 


1 


29 

A 
T 





_29 

1 


5 


1 


6 





36 


36 


36 





216 





216 





216 





1296 





1296 


6 


1 





6 


36 


36 


-36 








216 





216 





-216 


1296 





-1296 


7 


1 


-6 





36 


36 


36 





-216 





-216 





216 





1296 





1296 


8 


1 





-6 


36 


36 


-36 








-216 





-216 





216 


1296 





-1296 


9 


1 


V2 


V2 


4 


4 





2 


Ay/2 


Ay/2 


Ay/2 


Ay/2 








16 


8 





10 


1 


-V2 


V2 


4 


4 





-2 


-Ay/2 


4^/2 


-A-y/2 


Ay/2 








16 


-8 





11 


1 


-V2 


-V2 


4 


4 





2 


-Ay/2 


-Ay/2 


-Ay/2 


-Ay/2 








16 


8 





12 


1 


V2 


-V2 


4 


4 





-2 


Ay/2 


-Ay/2 


Ay/2 


-Ay/2 








16 


-8 
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3 

v/2 


3 


9 


9 





9 

2 


27 
V2 


27 
V2 


27 
V2 


27 
sft 








81 


81 
2 
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1 


3 

V2 


3 

V2 


9 


9 





9 
2 


27 
V2 


27 
V2 


27 
V2 


27 
V2 








81 


84 
2 
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1 


3 

V2 


3 


9 


9 





9 
2 


27 


27 


27 


27 








81 


81 
2 
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1 


3 

y/2 


3 

V2 


9 


9 





9 
2 


27 
y/2 


27 


27 
V2 


27 
V2 








81 


84 
2 






f eg = (pJ*Jy,e>,fc,f?,--.,f3) T , where ft* = 2P + (j 2 + f y )/p, ft* = (.,;. j 2 ) ,>. 
f? = 3*3y/P, fs q = V + 2P)j x /p, ft q = (e> + 2P)j y /p, fa = (4P + fjp + j 2 /p)j x /p, 
AT = (±P + j 2 JP + j 2 y /p)jy/P, fi! = (2P + 3 2 Jp-3 2 y /p)3z/p, = (-2P + J 2 jp-J 2 /P)3y/P, 
fH = 2{b+2)pR^HM)RTUl+3 2 y)/P+Ul+3 2 y ) 2 /p", Al = i(b+4)P+(3 2 x +3 2 )/p]j x 3y/p 2 , 
fa = [(& + 4)P + (j 2 + j 2 y )/p](j 2 x - j 2 )/p 2 , where P = pRT, and e' = bpRT + fjp is the 
twice of total energy e. The recovered Navier-Stokes equations are as follows: 
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dj x | d 
dt dx 



d_ 

P J dy 
d pRT du y 



dP 

dx dy 1 s 7 



Jxjy 
P 

du. 



dx 



+ 7?)] 
dy 

d pRT 2 du x du y 
dx s 5 b dx dy 



pRT du x du y 



sq dx dy 



)]• 



(48b) 



9jy 

dt 



d_ 

dx 



Jxjy 
P 



, i 

p 



d_ 

dy 



dP + d ^pRT { du % 
dy dx s 7 



du. 



dx 



4[^d 

dy s 5 



dy 

2 , . du x du 



)] 



b dx dy 



pRT du x _ du y 
s & dx dy 



(48c) 



96 +^-[(e + P)j x /p] + -^[(e + P) h /p] 



dt dx 



dx s% 2 dx 
dy s 9 2 dy 



(2 



du x 2 du x 2 du y . 

K + ( 



dx 

(2 dUy 



b dx b dy 
2 ch^ 2 <9u. 



+ -l^)Uy\} 



dx 
,du. 



dy 

du x 
-I 

9?/ 6 <9x b dy ' y ' y dx dy 
When s 5 = Sq = s 7 = s 8 = Sq, the above Navier-Stokes equations reduce to 



u x }}. (48d) 



<9p ^Ta 

<9t <9x Q 



djg g (jejp/p) 

dt dxR 



= 0, 
dP 



p 



<9e 9 
dt dx a 



[(e + PR] 



_d_ 

dxe 



dx ~ a ^ 
dT 



2 + 



P a /3 u a 



(49a) 

(49b) 
(49c) 



where 



pRT 



p B = (2/3 - 2/6) 



pPT 



a/3 



du a dug 

*^7 + 



2 9m 



x A 



<9w 



5 Q/3 , (a,/3,7 = x,y). 



■;p dx a 3dx x j ^ x 
Similar to the case of the first MRT model, the second one works also well for shocked 



compressible fluids. For more general cases, similar modifications to the collision operators 
of the moments related to the energy flux should be made 45]. 
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IV. SIMULATIONS ON HYDRODYNAMIC INSTABILITIES 



Hydrodynamic instabilities are ubiquitous in natural and industrial processes. The 
Rayleigh- Taylor (RT) instability, Richtmyer-Meshkov(RM) instability and Kelvin-Helmholtz 
(KH) instability are highly concerned in weapon physics and inertial confinement fusion. For 
example, during the spherical implosion procedure, the high pressure applied at the outside 
of the shell drives a very strong shock wave towards the centre of the device. This shock wave 
first accelerates the interface to a high velocity. Towards the end of the implosion the inter- 
face is decelerated by a combination of shock waves reflected from the center of the device 
and continuous deceleration due to build up of high pressure in the thermonuclear material. 
Such a very complicated acceleration/deceleration behavior results in two processes, RT in- 
stability and RM instability. Since the implosion is generally not perfectly symmetrical, the 
shear at the interface induces the third process, KH instability. Hydrodynamic instabilities 
in such procedure influence significantly the implosion physics and weapon performance. 
In this section, we summarize our recent attempts on LB simulations on KH 5l| and RM 



instabilities 



42 



431 ] . When studying the RM instability, the system must be compressible. 



In the case of KH instability, the system can be compressible or nearly incompressible. As 
a first step, we attempted the case with nearly incompressible fluids. 



A. Richtmyer-Meshkov instability 

The RM instability arises when a shock wave interacts with an interface separating two 
different fluids. It combines various compressible phenomena, such as shock interaction 
and refraction, with hydrodynamic instability, including nonlinear growth and subsequent 
transition to turbulence, across a wide range of Mach numbers. In inertial confinement 
fusion, the RM instability causes mixing between the capsule material and the fuel within, 
limiting final compression and thus the ability to achieve energy break-even or production. 
The RM instability problems in the plane occur when a shock wave travels from a light 
medium to a heavy one or when the shock wave travels from a heavy medium to a light one. 
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1. Shock wave from light to heavy media 



A practical example for this case is that the shock wave travels from air to SFq. For such 
a case, in our LB simulations we set the following initial physical field, 

' (p, u x , u y ,p) l = (1.34161, 0.361538, 0, 1.51332), 

(p,U x ,U y ,p) m = (1,0,0, 1), 

(p,u x ,u y ,p) r = (5.04,0,0, 1), 

where the subscripts /, m, r indicate the left, middle, right regions of the whole domain. 
Such an initial configuration can be explained as below: the interface of the middle and 
right regions separates the light and heavy media; the interface of the left and middle 
regions is the shock front; the shocked light medium is in the left and the pre-shocked is in 
the middle regions. Initially, the two media have the same pressure and different densities 
and temperatures. The corresponding Mach number of the shock wave traveling from left 
is 1.2. The shock wave will hit the interface with an initial sinusoidal perturbation. The 
initial sinusoidal perturbation at the interface reads x = 0.25 x N x x dx + 0.008 x cos(207n/), 
where the cycle in ^-direction of initial perturbation is 0.1, the amplitude is 0.008, N x is grid 
number, and dx is grid size. The following boundary conditions are imposed: (i) inflow at 
the left side; (ii) outflow at the right side, and (3) periodic in the ^/-directions. 7 = 1.4 in 
the whole domain. 

Since the Mach number is 1.2, the compressibility effects in this case is not negligible. 
Figure 7 shows the density and pressure contours at four different times, t — 0, 0.06, 0.3 and 
1.15. When the shock wave passes the interface from the left, a reflected shock wave to the left 
and a transmission wave to the right are generated (clearly seen in pressure field at t = 0.06). 
The transmission wave has a certain curvature at this time. Due to the compression, the 
interface produces a small deformation, and the perturbation amplitude reduces slightly. 
At t = 0.3, the reflected shock wave has been out of the computational domain, and the 



transmission wave becomes flat, which is consistent with the theoretical analysis of [52 ] . 
The perturbation amplitude begins to increase under the pressure gradient, producing the 
bubble and spike structures. The misalignment of pressure and density gradients causes 
a deposition of vorticity at the top of spike structure, and a mushroom shape is formed 
eventually. Fig.8 shows the changes of perturbation amplitude and growth rate with time. 
The amplitude is defined as half of the maximum distance between the crest and trough. 
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FIG. 7: Snapshots of RM instability (from light to heavy medium): density and pressure contours 
at t = 0, t = 0.06, t = 0.3, t = 1.15, respectively. From deep to light color, the level corresponds 
to the increase of values. 



From Fig. 8 one can clearly find the initial decrease of perturbation amplitude. During this 
initial period, the growth rate is negative. 

Now we go to some theories to explain and validate the simulation results. Richtmyer 



53| proposed an impulsive model in the case of a reflected shock wave via modifying the 
linear theory of Taylor for Rayleigh- Taylor instability. According to the impulse model, the 
growth rate reads, 

da , A 

ai 



da ... 
— = kAuAidi, 
dt 



«o(l - -p-, 



where k = 2tt/X is the wave number, Am is the velocity change across the interface, Ai 
is the post-shock Atwood number, a\ represents the post-shock amplitude, ao is the initial 
amplitude, D denotes the incident shock speed. Cmpr = 1 — Au/D is compression ratio. 
According to the initial conditions, the solution is Cmpr = 0.84, da/dt = 0.063. In the 
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FIG. 8: Amplitude and growth rate changes with time (from light to heavy medium). 



experiments of Meshkov 54j and Benjamin 55|, the measured growth rates are only about 



one half of that predicted by the impulsive model. Zhang and Sohn [56[ developed a model 
for the growth of RM unstable interface from early to late times in the case of light-heavy 
transition. The amplitude growth reads 

da v 

~dl ~ 1 + k 2 v a x t + max[0, (ka^ 2 - (A^ 2 + 0.5}(kv t) 2 

where vq = kAuA\a\. As shown in Fig.8, the LB result for growth rate qualitatively agrees 
well with that of Zhang-Sohn model. The amplitude reaches the minimum value 0.0065 at 
time t = 0.05, so the compression ratio obtained in simulation is Cmpr = 0.0065/0.008 = 
0.81. By the least squares fitting, the growth rate of amplitude 0.03 is obtained, which is 
about one half of the growth rate predicted by the impulsive model and consequently is in 
good agreement with the experimental result. In the nonlinear stage, the simulation results 
agree qualitatively well with the perturbation model proposed by Zhang and Sohn. 



2. Shock wave from heavy to light media 

A practical example is that the shock wave travels from air to He. To better understand 
such a case, in our LB simulation, we set a planar shock wave with the Mach number 2.5 
impinging on a sinusoidal perturbation x = 0.1 + 0.008 x cos(207r?/), where the cycle and 
amplitude of initial perturbation are the same with the case where shock wave travels from 
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light to heavy media. The initial physical field is as below: 

r 

(p,u x ,u yj p) l = (3.33333,2.07063,0,7.125), 
< (p,u x ,u y ,p) m = (1,0,0, 1), 

(p,u x ,u y ,p) r = (0.138,0,0, 1), 

v 

The boundary conditions in the y-direction and at the left side are the same as the case where 
shock wave travels from light to heavy media. Two different boundaries are applied at the 
right side: outflow condition (case I) and reflecting boundary (case II). The computational 
domain is a rectangle 0.6 x 0.1 for case I and 0.3 x 0.1 for case II, respectively. 

Figure 9 shows the simulation results for density field. Figure (a) corresponds to the 
outflow boundary and figure (b) corresponds to the reflecting boundary. Here 7 = 1.4. The 
collision parameters in case I are S5 = 10 4 , 10 5 for the others, and in case II are S5 = 10 3 , 10 5 
for the others. Simulation results show the following physical procedure: When the shock 
wave passes the interface, a reflected rarefaction wave to the left and a transmission wave 
to the right are generated. The pressure of heavy fluid near the crest is greater than the 
light fluid pressure. Driven by the pressure gradient, the perturbation amplitude decreases 
with the interface motion to the right. Then, the peak and valley of initial interface invert, 
the heavy and light fluids gradually penetrate into each other as time goes on, the light 
fluid "rises" to form a bubble and the heavy fluid "falls" to generate a spike. In case I, the 
transmission wave continues to move to the right, and no longer interacts with the interface. 
The disturbance of the interface continues to grow, eventually forming a mushroom shape. 
In case II, the transmission wave reaches the solid wall on the right and reflects to the 
left, encounters the interface again. This is known as the "reshocking" process. Following 
reshocking, the interface is compressed, as seen from the kink in the bubble. Furthermore, 
the amplitude grows more rapidly than prior to reshocking, the increased growth is due to 
the additional vorticity deposited on the evolving interface during reshocking. The pressure 
contours and velocity vectors near the interface at time t = 0.08 are shown in Fig. 10. Figure 
11 shows the change of disturbance amplitudes with time, corresponding to case I and case 
II, respectively. Because of the reshocking effect, a significant difference between Fig. 11 (a) 
and Fig. 11(b) can be observed. 

The interface reversal phenomenon is observed in the second case. With the interaction 
between shock wave and interface, disturbance will grow continuously. In the early stage, 
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FIG. 9: Snapshots of RM instability (from heavy to light medium), (a) Outflow boundary. From 
top to bottom, t = 0, 0.02, 0.08, 0.16, respectively, (b) Reflecting boundary. From top to bottom, 
t = 0, 0.02, 0.04, 0.08, respectively. From deep to light color, the level corresponds to the increase 
of density. 




FIG. 10: Pressure contours and velocity vectors at time t = 0.08 (from heavy to light medium, 
reflecting boundary). From deep to light color, the level corresponds to the increase of pressure. 

logarithm of growth rate is nearly linear with time, while changes into the non-linear in the 
late stage, spikes and bubbles occur. 

B. Kelvin-Helmholtz instability 

During the later stage KH instability strengthens the nonlinear developing of RT and RM 
instabilities, enhances the small scale mixing. In some cases, it may break the spkies. But in 
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FIG. 11: Amplitude change with time (from heavy to light medium), (a) Outflow boundary, (b) 
Reflecting boundary. 

some cases, we failed to observe the full effects of KH instability. For example, in the Eagle 
Nebula, why has the famous "Pillars of Creation" so large scale structures, instead of being 
broken by many small scale vortices? There must be some mechanisms to restrain the KH 
instability. Therefore, people study the KH instability from two sides. How does the KH 
instability evolve? How to enhance or restrain the KH instability? The strong nonlinearity 
and multiscale interactions make difficult theoretical study. The very complex 3D behavior 
challenge experimental diagnosis. Our LB modeling and simulation aim to help understand 
better the KH instability from both the two sides. 

For investigating the Kelvin-Helmholtz instability, we set the following initial physical 
field, 

rf l) = /!L^_2Ll£5t«4(Jl), (50) 

v(z) = VL + 2 VR - VL 2 VR tanh(-^-), (51) 

P L = P R = P, (52) 

where we have two characteristic length scales, D p and D v , which are the widths of density 
and velocity transition layers, respectively, pi = 5.0 (pn = 2.0) is the density away from 
the interface of the left (right) fluid. V£ = 0.5 (vr = —0.5) is the velocity away from the 
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interface in ^/-direction of the left (right) fluid, and Pl {Pr)= 2.5 is the pressure in the 
left (right) side. The system can be approximately thought of as "incompressible". The 
whole calculation domain is a rectangle with length 0.6 and height 0.2, which is divided into 
600 x 200 uniform meshes. A simple velocity perturbation in the re-direction is introduced 
to trigger the KH rollup and it is in the following form 

u = uq sin(%) exp(— kx), (53) 

where Uq = 0.02 is the amplitude of the perturbation. Here, k is the wave number of the 
initial perturbation, and is set to be 107r. The time step is At = 10~ 5 . 

At the initial linear increasing stage of KH INSTABILITY, the amplitude t] of perturba- 
tion evolves according to the following relation, rj = r/oe 7 ', where 7 is the growth coefficient 
and is dependent on the gradient of tangential velocity and gradient of density around the 
interface. In other words, 7 is dependent on the width of velocity transition layer D v and 
width of density transition layer D p . We discuss separately the KHI in three cases, (i) D v 
is variable and D p is fixed, (ii) D p is variable and D v is fixed, (iii) both D p and D v are 
variable. The increasing rate 7 for cases (i), (ii) and (iii) are referred to as j v , 7 P and 7^, 
respectively. We numerically obtain j v , 7 P and 7^ via fitting the curves of In E x \ max (t) ver- 
sus the time t, where E x \ max (t) is the maximum of E x (x,y,t) in the whole computational 
domain, E x (x, y, t) = p(x, y, t)u 2 (x, y,t)/2 is the perturbed kinetic energy at the position (x, 
y) at each time step t. 

Although viscosity damps the evolution of the KH INSTABILITY, here we focus on 
cases such as in inertial confined fusion where effects of the viscosity are generally negligible. 
Therefore, throughout the simulations, r is set to be 10~ 5 to reduce the physical viscosity. 
Boundary conditions are as below. Periodic in the ^/-direction and outflow (zero gradient) 
in the x-direction. 

1. Velocity gradient effect 

Figure 12 shows the evolution of the density field for the case with D v = 4 and D p = 8 at 
four different times. At t — 0.3 the interface has been wiggling under the initial perturbation 
and velocity shear. A nicely rolled vortex occurs and develops around the initial interface 
after the initial linear growth stage. The vortex becomes larger with time and a mixing layer 
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FIG. 12: (Color online) Density evolutions of KH INSTABILITY simulated using the LB model, 
where D v = 4 and D p = 8,t = 0.1 in (a), t = 0.3 in (b), t = 0.5 in (c), and t = 0.7 in (d). 

forms around the initial interface. 

To investigate the velocity gradient effect, we fix the width of the density transition layer. 
Figure 13 shows the density field for various D v at the same time, where D p = 8, t = 0.6 
and D v = 4, 8, 12, 16 in (a)-(d), respectively. Five contour lines are plotted in each plot. It 
is clear that the width of the velocity transition layer significantly affects the evolution of 
KH instability. The larger the value of D v , the weaker the KH instability, and the later the 
vortex appears. In Figs, (a) and (b), large vortices have been formed demonstrating that the 
evolution is embarking on the nonlinear stage. While in Figs, (c) and (d), the evolution is 
in the weakly nonlinear stage. Figures (a)-(d) show that a wider velocity transition zone is 
helpful for stabilizing the KH instability. 

The peak kinetic energy E x \ max partly indicates the interacting strength of two different 
fluids. Figure 14(a) shows that logarithm of E x \ max versus time. The initial state shows 
a linear behavior. The slope k increases with decreasing the width D v . After the initial 
stage, ln(E x \ max ) approaches a saturation value via a nonlinear growth stage. During the 
initial linear stage, we have E x oc u 2 oc (e 7 *) 2 . So, the slope k here can be used to calculate 
the growth coefficient 7 in the linear growth stage, k = 27. The logarithm of 7 decreases 
linearly with D v [see Fig. 14(b)]. Our LB results confirm the theoretical analysis of Wang, et 
al. (57J. In the classical case, the linear growth rate is 7 C = ky/pip 2 (vi — V2)/(pi + P2) oc Av, 
where Av is the shear velocity difference. A wider transition layer decreases the local or 
the effective shear velocity difference Av, which results in a smaller linear growth rate and 
a longer linear growth time. 
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FIG. 13: (Color online) Vortices in the mixing layer as a function of D v at t = 0.6, where D v = 4 
in (a), D v = 8 in (b), D v = 12 in (c), and D v = 16 in (d). The density transition layer D p is fixed 
to be 8. 
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FIG. 14: (Color online) (a) Time evolution of the perturbed peak kinetic energy ^ x | max along 
the x-axis in ln-linear scale for various widths of velocity transition layer. The dash-dotted lines 
represent the linear fits to the initial linear growth regimes, (b) Linear growth rate as a function 
of the width D v of velocity transition layer. 
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FIG. 15: (Color online) (a) Time evolution of the logarithm of the peak kinetic energy E x \ max along 
the x-axis for various widths of density transition layers, (b) Linear growth rate as a function of 
the width D p of density transition layer. 

2. Density gradient effect 

The density gradient effect is investigated in a similar way. Here D v is fixed. The 
initial conditions are described as, (/?£, Vl, Pl) = (5.0, 0.5, 1.5) and (pr, Vr, Pr) = (1.25, 
—0.5, 1.5). Figure 15(a) shows evolution of the logarithm of peak kinetic energy E x \ max 
along the re- axis versus time t for various widths of density transition layers. Here D v = 2, 
Ax = Ay = 0.002, At = 10^ 5 . Results for D p = 0, 2, 4, 6, 8, 10, and 12 are shown. 
For fixed width of velocity transition layer and density difference, the linear growth rate 
first increases with the width D p . But when D p is large than a critical value which is 
about 6, it does not vary significantly any more [see Fig. 15(b) ]. During the linear growth 
stage, 7 P increases linearly with the logarithm of D p . Figures 14 and 15 indicate the effective 
interaction width of D p is less than that of D v . The LB results here confirm also the 
theoretical analysis of Wang, et al. j^]. In the classical case, the square of the linear growth 
rate is >y 2 = k 2 p 1 p 2 (v 1 - v 2 ) 2 /(pi + p 2 ) 2 oc (1 - A 2 )Av 2 , where A = (pi - p 2 )/(pi + P2) is 
the Atwood number. A wider density transition zone reduces the Atwood number around 
the interface. Then in the process of exchanging momentum in the direction normal to 
the interface, the perturbation can obtain more energy from the shear kinetic energy than 
in cases with sharper interfaces. Therefore, a thinner density transition layer is helpful to 
restrain the KH instability. 



33 



5 10 15 20 

D 

P 



FIG. 16: (Color online) The linear growth rate versus the width of density transition layer for 
R = 0.5, 1, 2, and 5. The initial density, shear velocity and pressure of the two fluids are (pi, vl, 
P L ) = (5.0, 0.5, 1.5) and (p R , v R , P R ) = (1.25, -0.5, 1.5). 

3. Hybrid effects of velocity and density gradients 

In practical systems, at the interface of two fluids with a tangential velocity difference, 
both the velocity and the density gradients exist. There is a competition between effects of 
the two kinds of gradients. We introduce a coefficient R = D p /D v through which we analyze 
the combined effects. The linear growth rate versus D p under various values of R is shown 
in Fig. 16. Here R = 0.5, 1, 2, and 5, as shown in the legend. On the whole, the hybrid 
effect of the two kinds of gradients is to reduce the linear growth rate 7^. Only at small D p 
and when R > 1, the hybrid effect makes larger the linear growth rate. This indicates again 
that the effective interaction width of the velocity transition layer is wider than that of 
density transition layer D p . 
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V. CONCLUSIONS 



Both the LB and the hydrodynamic equations are simplified dynamic models of practi- 
cal systems. Compared with the latter, the former puts the physical modeling on a more 
fundamental level. When numerically study a physical procedure, the working dynamic 
model is not the one evolving continuously in space and time but the one discretized in 
the code. Improving the discrete template and reasonably adding viscosity term are in fact 
some remedies to the working dynamic model. Compared with the LB based on BGK ap- 
proximation, the MRT-LB introduces a new framework where various physical modes can 
be considered separately. The developed SRT-LB and MRT-LB are complementary from 
the sides of convenience and applicability. Compared with the hydrodynamic descriptions, 
both the SRT-LB and MRT-LB present new measurements for the deviations of systems 
from their thermodynamic equilibria. The LB model is being extended to study the com- 
pressibility effects, effects of shocking and detonation, thermal effects on the hydrodynamic 



instabilities 



51] and multiphase flows 
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60(, etc., which are all-important issues in science 



and engineering. 
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